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ABSTRACT 

A considerable number of astrometric binaries whose positions on the sky do 
not obey the standard model of mean position, parallax and linear proper motion, 
were observed by the Hipparcos satellite. Some of them remain non-discovered, 
and their observational data have not been properly processed with the more 
adequate astrometric model that includes nonlinear orbital motion. We develop 
an automated algorithm based on "genetic optimization", to solve the orbital 
fitting problem in the most difficult setup, when no prior information about the 
orbital elements is available (from, e.g., spectroscopic data or radial velocity 
monitoring). We also offer a technique to accurately compute the probability 
that an orbital fit is bogus, that is, that an orbital solution is obtained for a 
single star, and to estimate the probability distributions for the fitting orbital 
parameters. We test this method on Hipparcos stars with known orbital solutions 
in the catalog, and further apply it to 1561 stars with stochastic solutions, which 
may be unresolved binaries. At a confidence level of 99%, orbital fits are obtained 
for 65 stars, most of which have not been known as binary. It is found that reliable 
astrometric fits can be obtained even if the period is somewhat longer than the 
time span of the Hipparcos mission, i.e., if the orbit is not closed. A few of the new 
probable binaries with A- type primaries with periods 444-2015 d are chemically 
peculiar stars, including Ap and A Boo type. The anomalous spectra of these 
stars are explained as admixture of the light from the unresolved, sufficiently 
bright and massive companions. We estimate the apparent orbits of four stars 
which have been identified as members of the ~ 300 Myr-old UMa kinematic 
group. Another four new nearby binaries may include low-mass M-type or brow 
dwarf companions. Follow-up spectroscopic observations in conjunction with 
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more accurate inclination estimates will lead to better estimates of the secondary 
mass. Similar astrometric models and algorithms can be used for binary stars 
and planet hosts observed by SIM PlanetQuest and Gaia. 

Subject headings: astrometry — binaries: general 

1. Introduction 

One of the most promising and presently available ways to discover invisible compan- 
ions (including brown dwarfs and giant planets) to nearby Galactic stars is to analyze high- 
accuracy astrometric positions of the latter for the presence of the reflex Keplerian motion 
caused by an orbiting companion. In the relatively near future, with the advance of such 
space-borne astrometric instruments as SIM and GAIA, it will become one the main instru- 
ments in the search for habitable smaller planets, too. For the time being, the capabihties of 
the method are limited by the moderate precision of the available astrometric data. The Hip- 
parcos Intermediate Astrometry Data (HIAD), at a single point precision of roughly 10-15 
mas, is just good enough to detect brown dwarfs around nearby solar-type stars, reliably so 
in conjuction with spectroscopic measurements (Halbwachs et al. 2000; Pourbaix & Jorissen 
2000; Torres 2006). Perhaps, giant planets of intermediate orbital periods could also be 
detected around a few nearest stars. 

Obtaining a robust orbital fit for an astrometric binary becomes increasingly difficult 
when the orbital period P exceeds the time span of Hipparcos measurements, which is about 
3.2 years. If the period is 6 years, still more than half of the orbit is represented in the 
data, and astrometric analysis may provide an independent solution. When the period is 
longer than 9 years, the astrometric data alone can provide only ambiguous results, since 
the almost linear segment of the orbit is hard to distinguish from the regular proper motion. 
Generally, the astrometric orbital fit is a non-linear 12-parameter adjustment problem that 
includes five astrometric parameters and seven orbital parameters (Appendix A). In some 
cases, e.g., long-period orbits, or highly inclined orbits aligned with the line of sight, the 
parameters become so entangled in the non-linear condition equations, that the fit becomes 
ill-conditioned. Therefore, any orbital solution should be supported by analysis of confidence 
intervals for all fitting parameters. If the astrometric solution is sufficiently well-constrained, 
approximate standard deviations of fitting parameters can be derived in the near vicinity 
of the solution point by covariance analysis, a recipe of which is given in Appendix B. Ill- 
conditioned solutions will have nearly singular normal matrices that can not be inverted, 
or large variances for some of the parameters. In more complicated cases, some of which 
we investigate in this paper, the confidence intervals and probability distributions can be 
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calculated by extensive Monte-Carlo simulations. When the optimization problem is almost 
singular with respect to particular nonlinear parameters (most often the eccentricity) , direct 
mapping of the objective function on a sufficiently fine grid of parameters is carried out. 

Another crucial problem that arises in fitting Keplerian motion for stars previously not 
known as binaries, is estimation of the probability of a null hypothesis, that is, that the star 
is single, and the detected perturbation in the astromctric residuals is caused by a chance 
occurrence of random noise in the data, or other effects, unrelated to binarity. Wc offer a 
robust and straightforward (if somewhat computationally heavy) method of confidence level 
estimation of a binary detection. 

These techniques are tested on a set of 1561 stars in the Hipparcos catalog with the 
so-called stochastic solutions. These data often represent gross errors or failures in the 
data reduction. Due to the type of detector used in the main instrument of Hipparcos 
(a grid of long vertical slits and a photomultiplier behind it with a non-uniform response 
across the pointing field of regard), many of these failed solutions originate in the wrong 
assumption about the target positions, or in the lack of knowledge of the relative position 
and brightness of visual binary components. Successful attempts have been made to use 
more accurate information about stochastic stars from ground based and Tycho-2 data sets 
to reprocess the published Hipparcos transit data, resulting in accurate solutions for a few 
hundred stars (Falin & Mignard 1999; Fabricius & Makarov 2000). The remaining stars 
with stochastic solutions arc prime suspects for yet unknown binaries. Orbiting pairs with 
separations smaller than 0.1 arcscc were not resolved by Hipparcos; on the other hand, 
a nonlinear apparent motion of the photocenter as small as several milliarcseconds could 
break the regular 5-parameter solution. Long-period orbits could be detected as additional 
acceleration of proper motions or as discrepant proper motions from the short-term Hipparcos 
data and the long-term Tycho-2 definitions (Makarov & Kaplan 2005). The overlap between 
such accelerating astrometric binaries and the list of stochastic objects is substantial, but 
presumably a lot of binaries on shorter orbits are left out of this analysis, as they describe 
a strongly nonlinear track during the 3.2-3.5 time span of Hipparcos observations. 



2. Using HI AD to compute orbits 

2.1. Correlated input data 

The star abscissae data in the HIAD was derived by the two Hipparcos data reduction 
consortia, FAST and NDAC, almost independently, but from the same observational data 
(ESA 1997, Vol. 3). Typically, each great circle measurement produced a pair of abscissa 
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data points, one derived by FAST, and the other by NDAC. These pairs of measurements 
arc statistically correlated. The correlation coefficient was estimated and published in the 
HIAD record, when appropriate. However, some of the great circle results were processed 
(or accepted) by only one of the consortia, in which case the correlation is zero. Since 
any astro metric or orbital fit is a least-squares adjustment minimizing the on abscissa 
residuals, internal correlations must be properly taken into account. 

In this paper, we obtain separate orbital fits based on NDAC and FAST data. In our 
judgment, these data arc derived from the same observations with the same instrument, and 
they should be strongly correlated. The correlations given in the HIAD arc probably strongly 
underestimated, and using a weighted combination of these data may have an adverse effect 
on the confidence estimation. 



2.2. Multiple solutions and local minima 

The orbital motion model described in Appendix A contains 3 nonlinear and 9 linear 
parameters. After Pq, e and Tq are fixed by a nonlinear optimization method, the rest of 
parameters {A, B, F, G, and 5 astrometric parameters) can be found by solving a linear 
system of equations. 

Solving a multidimensional optimization problem is difficult. Deterministic methods 
such as the widely used Powell minimization on starting conditions can converge to different 
local minima. Two-dimensional optimization is sometimes marginally solvable by brute force, 
i.e., by walking over a grid of initial values and initiating optimization from each of these 
starting values. For 3 and more dimensions, however, this approach is time consuming and 
a better approach is welcome (Section 6). 

3. Genetic optimization algorithm 

We employ a method described in Storn & Price (1995), which is a generalization of the 
genetic algorithms for continuous functions optimization. The differential evolution (DE) 
algorithm is capable of handling nondifferentiable, nonlinear objective functions. It is one of 
the methods used in Mathematica NMinimize routine. Although our objective function 
is differentiable everywhere except several singular points (Appendix B) and, with respect 
to a few parameters, linear, we have chosen this algorithm for its enhanced ability to handle 
complex surfaces in the parameter space with multiple minima, as discussed in Section 4. 

In a population of potential solutions in a n-dimensional search space, a fixed number 
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of vectors are randomly initialized, then evolved over time to explore the search space and 
to locate the minima of the objective function. 

At each iteration, called a generation, new difference vectors are computed by subtract- 
ing two vectors selected randomly from the current population (mutation). The difference 
vector is added to a third vector, which may be the current optimal vector or a randomly 
selected one. The resulting candidate is mixed with the current best vector, imitating mixing 
of chromosomes in sexual reproduction, where each coordinate component corresponds to a 
chromosome. The resulting trial vector is accepted for the next generation if it yields a reduc- 
tion in the vahie of the objective function. This process, imitating biological evolution, is not 
guaranteed to find the best solution that corresponds to the global minimum of the objective 
function, just as evolution does not always select the best possible mutation. However, the 
probability of finding the global minimum, or a very close solution to the global minimum is 
increasing with the population size. Therefore, at the expense of extra computing time, we 
can find the 'nearly' best solution within a given model, localizing most of the local minima 
(alternative solutions) at the same time by assuming large numbers of generation vectors. 
This allows us to investigate the multitude of alternative solutions, rather than rely on a 
single possible solution that yields a sufficiently small x^. 

Existence of multiple solutions with roughly equally reduced is a warning that the 
pattern of noisy data is untenable for a reliable binary solution, or that some other non- 
binarity effect is responsible for the observed perturbation. 

4. Distributions, probabilities and confidence intervals 

Instead of the standard F-test, unsuitable for strongly nonlinear models, the following 
Monte-Carlo test is used to evaluate the robustness of our orbital solutions. We generate 
a set of observations of a star with coordinates, parallax and proper motions correspond- 
ing to the best fit 5-parameter model at the same observation times. Normally distributed 
random numbers with variances corresponding to the estimated (formal) measurements' er- 
rors are added to each recorded transit time. After that the realization is reduced for both 
5-parameter and 12-parameter models, resulting in xl^ and Xi2i ior i — 1,2, . . . , N, respec- 
tively. The resulting fits are compared with the original, unperturbed fits with corresponding 
X?2o and xio statistics. 

We estimate the fraction of trials for which the ratio of the 5-parameter fit to that of 
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the 12-parameter fit exceeds that ratio for the actual data: 

1 ^ 

i=l 

where Q{x) is a threshold function, equal to 1 for x > 0, and for x < 0. The confidence 
of rejecting a null hypothesis is 1 — p. This test is numerically equivalent to the F-test for 
a linear problem at sufficiently high N. 

If the probability test is pssed, extensive Monte-Carlo simulations are carried out to 
estimate the errors for each parameter. Each estimated transit time is perturbed by a 
normally distribTitcd random number with a standard deviation equal to the formal transit 
time error (given in the HIAD) and the reduction process is repeated. This produces a 
distribution of true transit times given observation, provided we have accurate transit time 
error estimates, and the errors are truly Gaussian. By reducing these data, we obtain the 
distribution of underlying parameters given observation. This method of parameter error 
estimation is known as parametric bootstrap (Hastie 2001). 

At least 1000 trials are performed for each object. From the collected Montc-Carlo data, 
we build histograms for each fitting parameter such as e, Oa, Tq, which arc subsequently used 
to estimate 99, 95 and 68 percent confidence intervals. The distribution of solutions is often 
significantly different from a Gaussian bell curve. With increasing number of measurements, 
and decreasing relative errors of parameters, the parameter likelihood function becomes con- 
fined to an area where linear approximation is quite adequate for error estimation (Appendix 
B). 

Figure 2 shows computed for two representative stars, as functions of period P in days 
and eccentricity e, while keeping Tq and the remaining 9 parameters constant. The spacing 
between contours is 50, the darkest color corresponding to the smallest x^- We find a 
well defined global minimum for HIP 1366, while for HIP 20087 we fail to find a consistent 
solution. Both have multiple minima (not all of them evident in this two dimensional contour 
plot). Because of the non-uniform time cadence of HIP 20087 observations and a smaller 
number of measurements (26 vs. 32 for HIP 1366), the map for this star is more complex. 
Another feature emerging from these plots is that the minimization favors high eccentricity 
solutions. This is not a feature of this minimization algorithm, because any x^-based global 
optimization algorithm is biased toward high eccentricity in poorly conditioned problems. 
As noted by D. Pourbaix (2005)^, the effective number of degrees of freedom goes up at 
high eccentricity, making it possible to find a fit with anomalously small residuals if the 



^http://wwwhip.obspni.fr/gaia/dms/texts/DMS-DP-02.pdf/ 
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Fig. 1. — Period cumulative distribution for HIP1366. The total number of trials is 14700. 
Horizontal lines show bounds of 68% confidence interval. 
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number of observations is small. Possible solutions to this problem would be to introduce 
a penalty function for high eccentricity solutions or use Bayesian approach with a low prior 
probability of high eccentricity solutions. We decided against using these techniques in this 
work to keep the algorithm as simple as possible. Existing Hipparcos data do not preclude 
high eccentricity solutions; and rigorous estimation of prior eccentricity distributions is a 
project beyond the scope of this paper. 



HIP 1366 HIP 20087 
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Fig. 2. — contours for two stars as functions of period and eccentricity. 



5. Comparison with known binary solutions 

To verify the algorithm, we run it on a dataset of 235 stars with known orbital solu- 
tions in Hipparcos. The initial period, eccentricity and periastron time for these Hipparcos 
solutions were often adopted from spectroscopic data. In our work, we assume no prior 
knowledge of any orbit elements and compute the fits from scratch. The genetic algorithm 
is initialized with 30 points randomly distributed in phase space and differential evolution 
simulations are carried out for 100 generations for each star. The process is repeated 1000 
times for each star with input data modified by adding normally distributed random num- 
bers with standard deviations equal to the formal errors. The results of this verification 
analysis are shown in Fig. 3, where the period of orbiting binaries derived by us and their 
95% confidence intervals are put to comparison with the periods from the Hipparcos catalog. 
The genetic algorithm fails to converge for four long period stars HIP 32349, 5336, 37279 
and 95501. 

The majority of orbital solutions are in good agreement as far as the crucial parameter 
P is concerned. It is to a large degree an external check for our method, since the Hipparcos 
orbital solutions are largely based on more accurate spectroscopic orbits for the parameters in 
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Fig. 3. — Comparison of orbital periods estimated by the genetic evolution algorithm versus 
periods in the Hipparcos catalog of orbital solutions. Only FAST consortium data are shown. 
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common (most notably, for P and e) . It is recalled that no apriori information or constraints 
are used in our algorithm. With a few exceptions, the periods shorter than 1000 d appear 
to be reliable. Binaries with longer periods exceeding the observation time span (typically, 
3.2 yr) obtain uncertain solutions. Such systems will benefit from using spectroscopic or 
interferometric constraints. Four binaries with short periods in Hipparcos (P < 100 d) 
obtain much longer periods in our unconstrained fits. These could be hierarchical multiple 
systems, for which radial velocity monitoring tends to pick up the short-period, high velocity 
amplitude signal, while the astrometric method is more sensitive to long-period signals of 
large angular excursions. 



6. Comparison with traditional optimization methods 

To compare the performance of our genetic optimization algorithm with traditional 
optimization methods, we implement a brute force grid search optimization algorithm and 
apply it to the same set of 235 stars with known orbital fits in Hipparcos. The minimal value 
of is obtained by Powell optimization algorithm with initial starting points for eccentricity 
at and 0.5, 18 starting points for period spaced evenly in logarithmic scale between 0.1 
years and 30 years, and 4 starting points for Tq spaced evenly between and P. This 
number of starting points 18 • 4 ■ 2 = 144 corresponds to the lowest empirically determined 
density of the grid that still provides the highest possible rate of correct solutions. The 
Powell minimization is iterated until the relative change in is less than 0.01. The correct 
solutions are defined as solutions for which the fitted period is within a factor of 1.5 of 
the spectroscopically derived period. For the FAST consortium data, the grid optimization 
routine yields 144 correct solutions, requiring on average 19500 function estimations. 

The differential evolution optimization algorithm with 30 initial points and 30 gener- 
ations produces 152 correct solutions, with 140 solutions identical to the grid algorithm 
solutions. Each DE solution is followed by a Powell optimization on the remaining linear 
parameters with the same 0.01 tolerance on x^, totaling an average of 50 function estimates. 

In general, the DE global optimization algorithm requires approximately 20 times less 
function estimations to find the global optimum. Since the search space is only 3-dimensional, 
the advantage of the DE algorithm is not overwhelming. This algorithm, however, allows us 
to use CPU-intensive methods for estimating errors of free parameters, which is an important 
and subtle aspect of orbit optimization. It requires days of CPU time instead of months on 
one of authors laptop. For problems with larger number of dimensions, the DE algorithm 
may be even more attractive. It should be noted that the rather low density of starting 
parameter grid in the grid optimization technique stems from the moderate intrinsic accuracy 
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of Hipparcos data (a few milliarcseconds) , which is comparable to the detectable orbit size. 
The future space astrometry missions (SIM and Gaia) will operate at a factor of 10^ higher 
accuracies, dramatically increasing the density of the initial grid. 

7. Unconstrained solutions for stochastic stars 

We further indiscriminately applied the genetic optimization algorithm to all Hipparcos 
stars with stochastic solutions. Among the solutions demonstrating stable convergence, we 
selected 65 stars with confidence above 99%, estimated as described in section 4. A literature 
search for these 65 stars resulted in 11 systems with sufficiently accurate spectroscopic orbits, 
listed in Table 3. Astrometrically estimated periods arc usually quite close to the more 
accurate spectroscopic periods at P < 1000 d. A few discrepant cases arc present when the 
periods from FAST and NDAC data are different by a significant fraction of the true quantity. 
Generally, the FAST solutions seem to be more rehable. Probably, the FAST consortium 
used stricter criteria in treating statistical outliers, occasional observations deviating far 
from the true astrometric abscissa, than the NDAC consortium. Such outliers can disturb 
the nonlinear optimization process, generating multiple wrinkles and gradients in the 
function. But some of the discrepant cases find more interesting astrophysical explanation. 

The star HIP 84949 = V819 Her is a complex triple system. The inner binary is 
eclipsing of Algol type with a period of 2.23 d. The outer companion of G8IV-III type is 
an active spotty star rotating with a period of about 86 d (van Hamme et al. 1994). It 
can hardly be a coincidence that our best NDAC solution produced a period of 89 d. In a 
binary system with a variable component, the photoccnter moves along the line connecting 
the components synchronously with the light curve, the effect known as Variability Imposed 
Movers (VIM). The magnitude of the VIM astrometric excursion depends on the angular 
separation, magnitude difference and the amplitude of variability. Accidentally, the VIM 
effect in HIP 84949 generated an orbital fit with the smallest reduced with NDAC data, 
which was not at all a typical draw, since all the percentiles correspond to long-period fits, 
i.e., the real outer orbit. The FAST solutions are not affected by the optical variabihty. 
This case exemplifies the difficulties and hazards of purely astrometric orbital solutions for 
complicated multiple systems with variable components. 

Another piece of evidence that most of our discoveries are real long-period binaries is 
the high rate of occurrence of accelerating stars (Makarov & Kaplan 2005). This is not an 
independent check though, because apparent accelerations are perceived either from the same 
Hipparcos observations, or from a comparison of Hipparcos and Tycho-2 proper motions. We 
find 19 stars in common with the two catalogs of accelerating stars in (Makarov & Kaplan 
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2005). 

Finally, the extensive Geneva-Copenhagen (GC) spectroscopic survey (Nordstrom et 
al. 2004), and other radial velocity observations in the literature provide a truly external 
verification of the binary nature of some of our solutions. The GC survey only includes 
Hipparcos stars of F, G and K type. In some cases, only one RV measurement is available, 
so that orbital motion is impossible to establish. In most cases, only a few observations are 
available, which is not enough to derive a orbital fit, but is sufficient to detect orbital motion. 
Combining these data with Hipparcos astrometric data might be a subject for future work. 
Two indicators in the GC catalog arc especially useful in this respect: the probability of 
constant RV and a flag for spectroscopic binaries. Among the 54 new binaries from our 
analysis, 24 have these indicators of variable RV in the GC catalog. A few more stars have 
been known as spectroscopic binaries from earlier investigations, but without reliable orbit 
estimations. Such is the case with HIP 23221 = 63 Eri, which is orbited by a white dwarf 
companion, HIP 25732 = HD 3615 (Grenier et al. 1999), HIP 38146 = HD 63660, a 
com star (de Medeiros & Mayor 1999), HIP 43352 = HD 75605, an UMa member (de 
Medeiros et al 2002), HIP 62512 = HD 111456, another UMa member (Freire Ferrero et 
al. 2004), and HIP 76006 = HD 138525, an F6 giant (de Medeiros et al 2002). 



8. Overview of Astrometric Binaries 

The 54 previously unknown binary systems (Tables 1-2) range in spectral type from A2 
to MO. Some of the early type stars stand out by their chemical peculiarities, most notably 
the A-Boo type star HIP 32607 = a Pic, the Am star HIP 30342 = v Pic, and the Ap star 
HIP 1366 = HD 1280. 

• HIP 1366 = HD 1280. This star may also be quite young. Non-magnetic Ap stars 
are often found in SB2 systems, and their companions are usually magnetic Am stars. 
This makes the subtle discrimination of Ap and Am stars complicated for unresolved 
systems. Carrier et al. (2002) derived spectroscopic orbits for 16 Ap stars, four of 
which may be Am stars. They found orbital periods from 1.6-2422 d, which makes 
our solution, 1017 d, a rather common case. Generally, Carrier et al. (2002) find 
distributions of orbital elements similar to those of field main-sequence stars, except 
for the distinct deficit of systems with P < 3 d. The latter fact is explained in terms 
of the tidal synchronization in tight binaries, which helps to retain high rotational 
velocities preventing the occurrence of chemical peculiarity. From our solution based 
on FAST data, the visible semimajor axis of the photocenter is about 1 AU, implying 
a rather massive companion, perhaps also an A star. HD 1280 has not been identified 
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as SB. 

• HIP 30342 = HD 45229. Also known as Pic, this Am and chemically peculiar 
star has drawn considerable attention in the literature. Although Am stars are often 
found in SB2 systems (Abt & Levy 1985), ours is the first indication of the binary 
nature of this object. Using our FAST parameters for a, tt and P, we estimate the 
mass of the companion at only 1/4 of the primary mass. The FAST corrected mass 
function is (MI/M^^Jfast = 0.09. 

• HIP 32607 = HD 50241. Ahas a Pic, a well-studied, bright A Boo-typc star, not 
previously known to be binary. Our derived a is close to 1 AU, indicating a large 
mass ratio. This comphes with the suggestion of (Faraggiana et al. 2004; Faraggiana & 
Bonifacio 2005) that many (may be all) A Boo stars are unresolved binaries with com- 
posite spectra, which explain the observed chemical under abundances. The testable 
prediction here is that the unresolved companions have similar brightness to the pri- 
maries, lest their contribution to the observed spectrum is too small to account for 
the peculiarities of spectral lines. An alternative hypothesis of the A Boo phenomenon 
is based on the gas-dust separation in stellar envelopes (e.g., Andrievsky & Paunzen 
2003). 

A few primaries in our sample are nearby late-type stars, whose companions can be 
low-mass dwarfs from the bottom of the main sequence. A large number of M dwarfs and 
substellar objects (L and T type) are expected to roam the solar vicinity, but relatively 
few have been identified. The nearest low-mass stars, because of their intrinsic dimness, 
arc traditionally important objects of investigation. Substellar companions are of particular 
interest since the rate of such 'failed stars' in spectroscopic binaries is currently considered to 
be too low in comparison with the theoretical mass function expectation (the 'brown dwarf 
desert' problem). 

• HIP 38910. Our solutions for a and P from FAST and NDAC are discrepant for 
this K5V star, probably because of the very long period > 6 yr. It seems that only 
a small segment of the orbit is present in the Hipparcos data. The parallax is quite 
stable nevertheless, at 51-52 mas, that is 3-4 mas smaller than the value given in the 
catalog. Our NDAC-based solution seems the more reasonable for this system, from 
which a small mass ratio is estimated. The companion to HIP 38910 can be a nearby 
M dwarf. 

• HIP 84223 = GJ 1213. A nearby K7 dwarf with apparently a long period. The 
updated parallaxes are above 40 mas, which warrants keeping this star in the Gliese 
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catalog. With a period of several years and Ua probably less than 1 AU, the companion 
can be a late M dwarf. 

• HIP 107143 = GJ 836.3. This K3Vp dwarf has a stable solution with a„ = 33 mas 
and P 1350 d. The updated parallax stays just above 40 mas. The companion can 
be a late K or early M dwarf. 

• HIP 118212 = GJ 913. The latest star on our list, this nearby MO dwarf stands 
out by its updated parallax of 68 mas which is significantly larger than the cataloged 
value of 58 mas. Our orbital fit is relatively robust, suggesting a late M dwarf as a 
companion. 

Four early-type stars, HIP 12225 = HD 16555 = 77 Hor (A6V), HIP 32607 = HD 50241 
(A7IV), HIP 39903 = HD 68456 (F5V), and HIP 62512 = GJ 9417 (F5V) are prominent X- 
ray sources detected by Rosat. Such objects are occasionally found among the general field 
stars, as well as in nearby open clusters. Since the magnetic activity of the corona is believed 
to be suppressed in such massive stars because of the vanishing convective zone, it has been 
suggested that unresolved binary companions of smaller mass are the main contributors 
in the detected X-ray fiux. Finding more astrometric binaries bolsters this surmise. The 
latter star, HIP 62512 is of special interest. Beside the strong X-ray radiation, it is also 
an extreme UV source and a member of the ~ 300 Myr old UMa kinematic group. An 
astrometric excursion of roughly 1 AU coupled with a ~ 4 yr period we determine for this 
system imply a mass ratio of 0.5. Thus, the companion can be a young, hot white dwarf. 
Interestingly, we find three more UMa nucleus members in our sample, HIP 43352 = HD 
75605, HIP 49546 and HIP 61100 = GJ 1160, the latter is a known spectroscopic binary 
(section 5). With these new addition to the list of UMa binaries, the binarity rate in the 
nucleus of this group becomes uncommonly high. 

The star HIP 21965 = HD 30051 of spectral type F2/3 is suspected to be very young. 
Accurate age determination in the upper half of the main sequence is difficult. Our orbital 
fit suggests a mass ratio of about 0.3, so that the companion should have subsolar mass. If 
resolved, it can yield an accurate isochrone age for the system. 

For HIP 23221 = HD 32008 = 63 Eri we obtain a similar orbit to the previous 

system, with a slightly longer period P ^ 850 d. This popular star has in fact been known as 
a spectroscopic binary, but a poor orbit has been available in the literature. The secondary 
component is a white dwarf. Our orbital fit suggests a fairly small mass ratio. 
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9. Summciry 

We present a method for unconstrained optimization of double stars orbital parameters, 
based on an adaptation of the genetic evolution algorithm optimization and Monte-Carlo 
simulations for error estimation. This method can become a useful tool for estimating orbital 
parameters of binary stars and planet hosts for experiments like SIM, which are expected 
to observe a limited number of pre-selected targets. Experiments like Gaia are expected to 
observe much higher numbers of objects and requirements for computing time are likely to 
become more stringent. However, as the total processing time for the computations in this 
paper was about two weeks of CPU time on a medium range laptop, we expect the entire 
Gaia astrometric orbital estimation to become feasible on large scale scientific clusters in the 
near future. This "embarrassingly parallel" algorithm is ideally suited for running in a grid 
environment and it will take advantage of the ever increasing computer performance. We 
believe that simplicity and ease of implementation of our algorithm outweigh the relatively 
high CPU time requirements. 

The 65 stars from the Hipparcos stochastic solution annex arc good candidates for 
follow-up spectroscopic studies, which are expected to yield more accurate orbital parameters 
constraining the astrometric fits. 

FactUty: HIPPARCOS 



A. Lineeirized equations for orbital fits 

The Hipparcos Intermediate Astrometry Data (HIAD), the basic data set to derive astromet- 
ric orbits, contains partial derivatives of the star abscissa with respect the five astrometric 
parameters of the standard model in equatorial coordinates. 



di —dui/da^ (Al) 

d2 = dai/d5 (A2) 

4 = dai/dU (A3) 

(^4 = dtti/dfia* (A4) 

0^5 = dtti/dus (A5) 



where is the abscissa in ith observation of a given star, a and 5 are the equatorial coordi- 
nates, — a cos S, H is the parallax, and fia* — fJ-a cos S and /is are the orthogonal proper 
motion components. The abscissa is defined as a great circle arc connecting an arbitrary 
chosen reference zero-point and the star on a fixed great circle (close to the scan circle). 
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In the small-angle approximation, a linearized equation for the observed abscissa differ- 
ence Acj = Oobs ~ flcaic) can be written as 

L/Co L/C -J 

J •' 3 ■' 

where are the elements of the vector of seven orbital elements, e = [ca, P, e. To, 
Per usual, Oa is the semimajor axis of the apparent orbital motion (i.e., the motion of the 
photocenter) , P is the orbital period in years, e is the eccentricity, Tq is the periastron 
time, uo is the longitude of the ascending node in the plane of orbit in radians, Vt is the 
position angle of the node in the plane of projection, and i is the inclination of the orbit 
(zero for face-on orbits). In this equation, the notations a* and 5 were changed to x and 
respectively, to make them consistent with the traditionally used tangential coordinates for 
apparent orbits. Eq. A6 holds only in the vicinity of a certain point in the 12D parameter 
space {a, 5, 11, //q*, /z^, e}, as long as the corrections to these parameters remain small. 

The apparent motion of a binary in the plane of celestial projection is described by 
(Heintz 1978): 



X = A(cos E-e) + FVI - e2 sin E (A7) 
y = 5(cos E-e) + GVl - sin E 

where x and y are the tangential coordinates, and E is the eccentric anomaly related to the 
mean anomaly M by Kepler's equation 

M = 27r "^ ~ ^E- esin^;. (A8) 

The Thiele-Innes constants are related to the remaining orbital elements by 

A = aa{cosuj cos Q — sin uj sin Q cost) (A9) 
B = aa(cosc<Jsinr2 + sinc<Jcosr2cosi) 
F — aa(— sinci;cosf2 — cosa;sinOcosi) 
G — aa(— sin a; sin n -|- cos a; cos cos i), 

where is the apparent semimajor axis of the photocenter, uj is the periastron longitude, 
fl is the node and i is the orbit inclination (i = 90° is an edge-on orbit). 



B. Covariances with linearized equations 

Exphcitly, partial derivatives ^ and ^ in eq. A6 are 

fBll 



dei 

dx X 
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where 



dy_ 

dua 
dx 

dP 
dy 
dP 
dx 

dy 
de 
dx 
dTo 
dy 
dTo 
dx 

du 
dy_ 

du 
dx 

M 
dy 

dn 

dx 

'di 
dy_ 

di 



y_ 

dE 

[-A sin E + FVl - cos E] — 

dP 

, f)E 

\-B sin E + GVl - cos E] — 

dP 

..^ . ^dE. ^ sin E'fcos — e) 

-A(l + sm E—) + F , ^ '— 

^ de^ yr^^(l -ecosE) 

„, ^dE. ^ sinEfcos-E — e) 

-B(l + sin E—) + G , ^ — 

^ de' x/F^(l-ecos^) 

dE 

-A sin E + FVl - cos E\ — 

diQ 

dE 



BsmE + GVl - e2 cos - 

diQ 

dA /- dF 

— + VI - sm 

dou duj 



cos £^ — e 
cos E — e 
cos E — e 
cos £■ — e 
cos E — e 
cos £■ — e 



dA 

duj 
dB 

duj 



dB , dG 

-z h vl — e"' smii/— — 

dio duj 



dA 



dF 



- + v/r^sinE— 

dB , dG 

_ + VT^s.nE- 

dA , ^ dF 

-TT- + V 1 - 6"= sm E—- 

di di 



dB 



dG 



— + Vl - e2 sinE— , 



27r(T - To 



p2(ecosE- 1) 
sin£' 

1 — e cos £^ 

-27r 

P(l - ecos£;) 
F 



(B2) 



Qa sin u! sin Q, sin i 
—Qa sin cu cos Q sin i 
Oa COS LU sin O sin i 
—tta cos a; cos ^2 sin i . 

Using these relations, the design matrix D can be computed by eq. A6, and the co- 
variance matrix C — {D'^D)~^ is readily computed. Astrometric abscissae Oj for the same 
star may have different standard errors in HIAD, making it necessary to employ a weighted 
least-squares routine. 

The research described in this paper was in part carried out at the Jet Propulsion Lab- 
oratory, California Institute of Technology, under a contract with the National Aeronautics 
and Space Administration. This research has made use of the SIMBAD database, operated 
at CDS, Strasbourg, Prance. 



on 

OA 

'di 
dB 

du 
OF 

duj 
dG 
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Table 1. Orbital solutions from NDAC data. 
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Table 1 — Continued 
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Table 2. Orbital solutions from FAST data. 
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c:7+0.16 
"-"-0.13 


522_gs 


loolll'g 


15lllf« 


44+10 


nil 


38596 


27tf 


2566lgi82 





fiS + 0-27 
6S-0.20 


1298l««'^ 


+ -L 
2»»-17() 


99llf 


62ll« 


31+1 


38910 


79+84 


0407+3199 





09+0.50 
•^^-0.22 


2063l?oo« 


164ll«« 


1171- 


53l?l 


51+? 


38980 


15tf 


991+.\ll 





74+0.25 
'^-0.37 


439123^^ 


295111, 


io5iir 


76l| 


34+1 


39681 


2ltY 


1579li^ 





66+''-22 


103ll|?| 
12511^ 


6OI™ 


79+31 
'^-20 


72l« 


15122 


39903 


soil 


922III 





1 /+0.05 
-■-^-0.04 


33Q+16 
■'•-'^-334 


35l«2 


3211 


49+1 


40015 


7+18 
' -1 


245l| 
833+1']' 





55+0.43 
"^"-0.27 


I96I70 


1071^7 


1161^5 


11211'^ 


3311 


42916 


23+2 





Lfi riQ 

CQ-t-U.UO 

oy_o.o8 


384+^° 


1421^7 
208lll;i 


^ ^ .4-71 

104175 


150li, 


36+1 


43352 


9+1 


ii74i;?i 





tri+0.19 
■^^-0.14 


451+488 


195+156 


138+? 


13+1 


49546 


utl^ 







70+0.21 
'°-0.25 


466112,4 


328III4 


70+165 
'9-26 


84lg 


16+1 


50180 


35lf 


i838lii;;^ 





07+0.31 
•J' -0.18 


1041+*20 
-'-'^^-'--680 


1941^7., 


150+150 

J^'^^-ioi 


77+1 


30+1 


54424 


10+? 


112412™ 





55+0.18 
■^"^-0.17 


615lfo«o 


3OOIIL 


126+1?'^ 


12811J 


1511 


54746 


12+^' 


10181^0 





O.1+0.16 

"■^-o.ie 


166lif 


8511? 


127122^ 


70+9 


19+1 


61100 


31+1 


14201?°^ 





61+°12 
"J^-0.10 


201+?!^ 


176l?4 


5111^ 


571^ 


41+1 


62512 


45+1° 


1600l™;i 





9n+0.27 
^'J-0.13 


299li°f 


sil^^s 


16411° 


108+2 


39+1 


73440 




448+I3 





(55+0.34 
"°-0.26 


109+^1 


244+24 
^^^-163 


84+i?« 


6811.1 


30+1 


75401 


19+12 


9631^°^ 





47+0.39 
*'-0.29 


677+292 
'-347 


121+17 


132+g 


74l?, 


19+1 


76006 


s+_l' 


545111 





"="-0.33 


2201^5 


75lli« 


111+61 
^^J^-80 


571- 


Mil 


78970 


13+1 


87111^ 





97+0.53 
^' -0.16 


6631^8^ 


5.1 + 178 


100l|« 


61+™ 


22+? 


80884 


14+2 


10081^° 





4f5+0.17 
^0-0.13 


11 (.+710 

llU_gg 


29ll?76 


118lif 


64l? 


10+1 


84062 


11+? 


7381^1 





oQ+0.28 
•^^-0.19 


9QQ+366 


97+15 


8lti 


60l| 


3511 



-24- 



Table 2 — Continued 



84223 


30t^ 


3607l?g3 


0.731 


0.26 
0.32 


19581^470 


151+i?« 


looll^? 


110+1?^ 


40 


fi 
-1 


84949 




1604l|«5 


0.761 


0.17 
0.17 


451+^^^ 


i4oii;:^ 


62+1** 


60ir> 


15 


fi 
-1 


85852 


811 


903+1? 


0.30l 


0.23 
0.16 


6401S 


18111,^ 


111+99 
iii_7g 


1351- 


11 


fi 

-1 


88848 




2053l«?J 


0.731 


0.16 
0.20 


660lSil 


253+^ 
^^'^-178 


103+1^3 


78+1 
13011;] 


30 


fi 
-1 


90355 


11 + 2 


304+?j 


0.40l 


0.48 
0.19 


193+«» 


155lg 


7ll«i 


30 


V2 
-2 


93137 


13+^ 


874l^» 


O.61I 


0.29 
0.26 


549li«3 


28l^ 


loill.^ 


89+;; 


19 


fl 
~1 


94347 


13+1 


63411^ 


0.241 


0.13 
0.11 


345+?^ 


or, + 15b 


64i!i;; 


146+1" 


23 


fl 
-1 


94802 


18+1 


1037+27 


O.lll 


0.08 
0.06 




91+f« 


102l'^9 


130+1 


24 


fl 

-1 


97063 


16+2 


1376lti 


0.251 


0.32 
0.14 


848lts5 


48+lp 


96lg 


36+1" 


26 


fl 
-1 


97690 


91? 


1063l*>;^ 


0.231 


0.44 
0.13 


007+436 
■'0'-289 


47+170 

^'-16 


109l'^« 


60l« 


10 


f 1 
-1 


98375 


«-2 


284+1" 


O.77I 


0.22 
0.40 


LI <_ii7 


iqi + 16 
iji__^7,) 


0^ + 105 
80_27 


93l« 


20 


fl 
-1 


103287 


30+^ 


1347l?g 


O.35I 


0.38 
0.23 


O7O+406 
^'■^-661 
lZUi_7g4 


339_i78 


197+161 


88+? 


251^2 


104440 




2429l?;fS 


0.651 


0.11 
0.07 


17+180 


147l??,6 


91II 


51 


fl 

-1 


105969 




867l« 


O.42I 


0.30 
0.22 


0-14+386 
■^-■-^-137 


93llf 


iiolg 


130+1 


15 


f3 
-2 


107143 


35+1 


1346lgJ'' 


0.291 


0.23 
0.16 


1092l?,lf, 


276+1?., 


loeljf 


41+1" 


40 


f2 
-2 


110291 


oq+17 
•^^-26 


7941^!] 


O.99I 


0.00 
0.19 


562+1^^ 


2341-3 


87+1^'' 
llllt! 


99+1* 


15 


f3 
-1 


113638 


12+f 


3661" 


O.40I 


O.-'iO 
0.26 


186+77 


108l|i» 


6411^ 


31 


f7 
-8 


117622 


13l? 


lOTOl^Q^ 


O.44I 


0.23 
0.17 


o2f5+628 
•J^u_248 


33+1 


961^7 


119l« 


20 


fl 
-1 


118212 


38l| 


oQQ+438 


o.58l°;l« 


P-r,9+248 


112+i«8 


IO9I2I 


1131^ 


681^ 



Table 3. Comparison of astrometric fits with known spectroscopic orbits. 



NDAC 



FAST 



HIP 


c P [d] 




e 


12062 


1375 ± 227 





99 


19832 


678 ± 27 





23 


37606 


384 ± 11 





85 


38018 


544 ± 11 





48 


40015 


252 ± 7 





99 


61100 


1284 ± 118 





62 


73440 


435 ± 18 





50 


84949 


89 ± 845 





88 


85852 


916 ±61 





26 


88848 


4457 ± 776 





95 


105969 


887 ± 35 





23 



(M|/Af2 



P[d] 



(Mi /Ml,) 



Spectroscopic 



P [d] 



fiM) 



Ref. 





01 




55 


+ 


25 


-0 


17 




08 




34 




12 




11 




01 




56 


+0 


23 


-0 


15 




14 


±[! 


12 




09 




16 




30 




20 


t°o 


04 




21 




22 




18 



5.16 


973 ± 140 





64 


0.023 


668 ± 34 





19 


0.056 


386 ± 10 





85 


0.069 


555 ±8 





42 


0..54 


246 ±5 





38 


0.046 


1384 ± 227 





59 


0.0064 


450 ± 13 





44 


47.6 


1609 ± 462 





72 


0.069 


906 ± 46 





18 


0.25 


1460 ± 740 





43 


0.086 


859 ± 41 





32 



+0 


22 




28 




37 




18 


l°o 


10 




39 


+0 


12 


-0 


11 


+0 


44 


-0 


27 




12 


IE! 


09 


+0 


34 


-0 


16 


+0 


16 


-0 


16 


+0 


14 


-0 


14 


+0 


15 


-0 


20 


+0 


30 


-0 


21 



0.032 
0.023 
0.033 
0.068 
0.010 
0.033 
0.0064 
0.041 
0.056 
0.094 
0.12 



905 ± 12 
717 ±3 
380.6 ±0.1 
551.9 
243.8 
1284.4 
467.2 ±9.7 
2018.8 ±0.7 
903.8 ± 0.4 
2092.2 ±5.8 
878 



0.260 ±0.032 
0.074 ± 0.029 
0.73 ±0.012 
0.39 
0.42 
0.50 

0.2170 ±0.077 
0.672 ± 0.002 
0.072 ± 0.031 
0.765 ± 0.013 
0.13 



0.0262 ± 0.0038 



0.00078 ± 0.00018 

0.00351 ± 0.00031 
0.113 ±0.011 



(1) 
(2) 
(3) 
(4) 
(4) 
(3) 
(1) 
(5) 
(6) 
(7) 
(8) 



Note. — References; (1) Latham et al. (2002); (2) Halbwachs et al. (2000); (3) Setiawan et al. (2004); (4) Halbwachs (2003); (5) Scarfe et al. (1994); 
(6) Fekel et al. (1993); (7) Fekel et al. (2005); (8) Pourbaix & Jorissen (2000) 



